"""
蒙特卡洛模拟算法求解圆周率（Pi）

问题描述：
    使用蒙特卡洛模拟方法估计圆周率 \(\pi\) 的值。
    其基本思想是通过随机抽样来解决计算问题。
    
算法原理：
    将一个单位正方形内切一个半径为0.5的圆。
    正方形的面积是1，圆的面积是 \(\pi * r^2 = \pi * (0.5)^2 = \pi / 4\)。
    
    在正方形内随机生成大量的点，统计落在圆内的点的数量和总点数。
    根据几何概率，落在圆内的点的比例应近似等于圆面积与正方形面积之比：
    
    \(\frac{\text{落在圆内的点数}}{\text{总点数}} \approx \frac{\text{圆面积}}{\text{正方形面积}} = \frac{\pi/4}{1} = \frac{\pi}{4}\)
    
    因此，圆周率 \(\pi \approx 4 \times \frac{\text{落在圆内的点数}}{\text{总点数}}\)。

核心概念：
    - 随机抽样：通过大量随机样本来近似目标值。
    - 概率统计：基于事件发生的频率估计概率，进而推断未知量。

应用场景：
    蒙特卡洛模拟是一种强大的数值计算方法，广泛应用于：
    - 金融风险评估
    - 物理学中的粒子模拟
    - 工程可靠性分析
    - 统计推断与抽样
    - 优化问题（如MCMC方法）

作者：斯黄
日期：2025年3月4日
"""

import numpy as np
import matplotlib.pyplot as plt

# 蒙特卡洛模拟计算圆周率
def monte_carlo_pi(num_points):
    points_in_circle = 0
    
    # 用于可视化
    circle_points_x = []
    circle_points_y = []
    square_points_x = []
    square_points_y = []

    for _ in range(num_points):
        x = random.uniform(0, 1)  # 在[0, 1]范围内随机生成x坐标
        y = random.uniform(0, 1)  # 在[0, 1]范围内随机生成y坐标
        
        distance = x**2 + y**2
        
        if distance <= 1: # 如果点在半径为1的圆内
            points_in_circle += 1
            circle_points_x.append(x)
            circle_points_y.append(y)
        else:
            square_points_x.append(x)
            square_points_y.append(y)
            
    pi_estimate = 4 * points_in_circle / num_points
    return pi_estimate, (circle_points_x, circle_points_y), (square_points_x, square_points_y)

# 可视化函数
def visualize_monte_carlo(pi_estimate, circle_points, square_points):
    plt.figure(figsize=(8, 8))
    
    # 绘制正方形
    square = plt.Rectangle((0, 0), 1, 1, color='blue', fill=False, linewidth=2)
    plt.gca().add_patch(square)

    # 绘制1/4圆
    circle = plt.Circle((0, 0), 1, color='green', fill=False, linewidth=2)
    plt.gca().add_patch(circle)

    # 绘制落在圆内的点
    plt.scatter(circle_points[0], circle_points[1], color='red', s=1, label='Points in Circle')
    # 绘制落在正方形但不在圆内的点
    plt.scatter(square_points[0], square_points[1], color='gray', s=1, label='Points in Square (Out of Circle)')

    plt.xlim(0, 1)
    plt.ylim(0, 1)
    plt.gca().set_aspect('equal', adjustable='box')
    plt.title(f'Monte Carlo Estimation of Pi: {pi_estimate:.4f}')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.legend()
    plt.grid(True)
    plt.show()

if __name__ == "__main__":
    num_simulation_points = 100000
    pi_est, circle_pts, square_pts = monte_carlo_pi(num_simulation_points)

    print(f"使用 {num_simulation_points} 点估计的圆周率 Pi: {pi_est:.6f}")
    print(f"实际的圆周率 Pi: {np.pi:.6f}")

    visualize_monte_carlo(pi_est, circle_pts, square_pts) 